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Abstract 

Some statistical mechanics approaches to jammed granular media are 
based on ensembles in which the stress state of the system is externally 
set. This paper proposes a new phase space to describe the microstates of 
a granular packing compatible to a given macrostate. The nature of this 
phase space is analyzed, showing that the consideration of the allowed and 
forbidden regions of every configuration (i.e. geometrical pattern) could 
be a relevant factor for the establishment of its probability and, therefore, 
of the expected properties of the sample. 

Many combinations of forces acting on a particle can keep it in static 
equilibrium. Every set of forces could be considered equivalent to a mi- 
croscopic stress field, but the kind of interaction and the geometrical re- 
strictions mean that not all stress states can be represented by any set. 
Consequently every local configuration has its respective allowed and for- 
bidden regions in the phase space. As a result, some points of the phase 
space are degenerate, and the density of states strongly determines the 
expected statistical distribution in the thermodynamic equilibrium. It is 
shown how this function just depends on the deviatoric stress. A first 
analysis of two-dimensional (2D) arrangements is included to clarify this 
assertion. 

1 Introduction 

Edwards was the first to propose that a statistical mechanics approach might 
be feasible to describe dense granular media [1]. Assuming that granular sys- 
tems have entropy, it was claimed that the volume plays the role of energy 
(V-ensemble) [1, 2, 3, 4, 5]. There are other formalisms based on the energy of 
the whole system, such as [6, 7, 8]. Following this work, an approach based on 
considering the elastic potential energy of particles has recently been presented, 
referred to as the elastic energy approach [9] . This approach considers the stress 
state of the system. 

In this respect, in recent years, some researchers have suggested that the 
stress of the system should be considered together with the volume (V-F en- 
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semblc) [10, 11, 12, 13] or even alone (F-ensemble) [14, 15, 16, 12, 17, 18]. This 
can be done via the force-moment tensor, and this approach is referred to here 
as the force-moment approach. It uses the concept of angoricity. 

Both the elastic energy approach and the force-moment approach consider 
the external stress as a control parameter. In practice, the main difference be- 
tween them is whether the statistical weight of the points in the phase space 
decays as a linear or a quadratic function of the stresses, depending on the con- 
straints of the ensemble. However, this paper tries to show why it is not only the 
Boltzmann factor of the points of the phase space that determines the expected 
distribution in the thermodynamic equilibrium but also the degeneracy of the 
points. This is precisely due to the existence of different configurations with 
respective allowed and forbidden regions in the phase space. 



2 Revision of elastic energy and force-moment 
approaches 

The elastic energy approach [9] was an attempt of expressing the Edwards' first 
formalism [1] in terms of elastic potential energy, rather than volume. As a con- 
sequence, a theoretical framework based on stresses (or strains) and geometries 
was developed to explain the expected features of disordered arrangements of 
jammed granular media. 

The thcrmodynamical formalism was set up starting from the assumption 
that, in thermodynamic equilibrium, the expected macromechanical behavior of 
the sample always take the same value. In particular, if the equivalent stiffness 
of the system were supposed to be the same in all microstates compatible with 
a macrostate described by external pressure, the clastic potential energy also 
would be the same in all of them. The ergodic hypothesis would be guaranteed 
not by averaging samples in time but in cycles of stress and breakage, and conse- 
quently, those arrangements in mechanical equilibrium corresponding to global 
stiffness values different from those expected in thermodynamical equilibrium 
would be considered as states out of equilibrium. 

However, the physical restriction that the whole arrangement should be a solu- 
tion of the elastic problem was not rigorously imposed in [9] . This was done by 
first considering limit situations (Voigt and Reuss states) and then integrating 
over the allowed regions of every configuration and relating the arithmetic aver- 
age of the stress to the external pressure. Only crystal-like configurations were 
analyzed. 

The force-moment approach [17, 18] properly imposes that the assembly 
must fit a solution of the elastic problem and can be applied to arrangements 
with different values of equivalent stiffness and, therefore, of energy. This means 
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that, if particles are small in comparison with the scale of variation of the macro- 
scopic stress field , the latter must coincide with the average equivalent stress 
field over the whole system ov,-. is the stress field which fits a solution of 
the elastic problem with average equivalent values of constitutive relationships. 
In consequence it is possible to set up an ensemble based on the additive force- 
moment tensor = (JijV, with V being the volume of the system. Then, the 
Boltzmann factor used in the canonical ensemble depends on a linear function of 
the stresses, whereas the external control parameter, the tensorial temperature, 
is fixed by angoricity. 

The elastic energy approach showed that consideration of the different con- 
figurations and the allowed regions of the phase space can be useful to determine 
why some configurations are more likely than others and why some of them are 
barely present. The same methodology, revised and extended to general cases, 
is followed in this paper, since it is useful for a better definition of the random 
close packing [19] state or of the maximally random jammed state [20]: if the 
expected probability of each configuration in equilibrium were known, it would 
be possible to determine which the expected packing ratio of the sample is. 
However it is only possible if the constrains of the ensemble are clearly estab- 
lished in a suitable phase space. The objective of this paper is to show up how 
the nature of the phase space is. 



3 Equilibrium and phase space 

In this section it is proposed the r, T phase space (r, T have been chosen because 
the Latin word for stress is tensio) for those ensembles in which the stress state 
of the granular system is externally set. Since granular H-dimensional systems 
(H = 1,2,3) are characterized by both the strong interaction between particles 
and the inelasticity of collisions, the usual 2K./V-dimensional T phase space (de- 
scribed by the position and the momenta (pi, <?,) of the iV grains) together with 
the complex and dissipative dynamics make it difficult to apply the usual statis- 
tical mechanics frameworks. Therefore it is interesting to use a different phase 
space in order to compare packings and in order to apply statistical analysis. 

This is possible by taking advantage from the fact that sometimes granular 
arrangements must correspond to a solution of a macroscopic elastic problem. 
This is the case of static, jammed and compressed packings in which the size 
of the particles is small in comparison with the scale of variation of the macro- 
scopic stress field (for instance, this approach would not be appropriate for 
diluted granular gases or systems with few particles). 

In this way, one possibility is the analysis of the ensemble according to the 
interaction forces of the particles. It has been done within the framework of the 
force network ensemble theory [14, 15, 16, 12]. Jammed packings are usually 
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hyperstatic (i.e., the amount of force components is substantially larger than 
the number of force balance constrains) and they can be compared as a col- 
lection of coordinates in a ziV/2-dimcnsional phase space, being z the average 
coordination number. However it was also proved that force balance constraints 
on every particle reduce the dimension of the space. 

Other possibility is the use of the equivalent stress (or force-moment) fields 
associated to every particle as coordinates of the phase space, instead of forces: 
every arrangement in mechanical equilibrium with an external stress field can 
be compared as a collection of points in a stress r (or force-moment T) phase 
space, irrespective of where the particles are (in the real space) and which their 
respective neighbors are. Then it is possible to obtain the most probable dis- 
tribution of points which satisfies the restrictions of the ensemble (total energy 
or average equivalent force- moment field). It is being referred to this space as 
"phase space" in the sense of that it is the space of all possible states of these 
system but it is very important to emphasize that the laws of classical mechanics 
are not applicable there. 

Before presented in depth which the nature of the r, T phase space is, some 
points on the different Boltzmann and Gibbs approaches on the definition of the 
equilibrium must be clarified [21]. The Boltzmann approach needs to introduce 
some macroscopic variables (e.g. kinetic energy for ideal gases) and equilib- 
rium is defined by reference to the macrostate, which is uniquely given by a 
set of values of the macrovariables. Usual restrictions for ideal gases are those 
given by the microcanonical (N,V,E), canonical (N,V,T) or grand canonical 
(fi, V,T) ensembles 1 . The expected distribution in the equilibrium is obtained 
by maximizing the Boltzmann's entropy S cx In f2 2 , and it gives which the most 
probable distribution of an ensemble subjected to the same constraints is. The 
problem is not mechanical but statistical one, because the role of the dynamics 
is trivial for this definition: it is being compared, among the solutions of a par- 
ticular physical problem, the class that it is compatible with the macroscopic 
knowledge of a system. 

On the other hand, the starting point for the Gibbs approach is to suppose 
that the points in the phase space are distributed according to a probability 
density which is invariant under the flow of dynamics, i.e. this distribution 
function is a solution of Liouville's equation. Equilibrium is defined for a situ- 
ation where the probability density function is not an explicit function of time 
and it becomes a function of the global constants of motion. The statistical 
mechanical analogues of thermodynamic quantities are either fixed external pa- 
rameters, related to phase functions or functionals of the density function (e.g. 
the Gibbs entropy S oc p In p) . 



1 V is the volume, E is the energy, T is the temperature and fi is the chemical potential. 
2 fl is the number of microstates compatible with a macrostate of the system 
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In the case of ideal gases, both approaches are successful and they give the 
same expected distribution for the same constraints. Moreover, in these thermal 
systems particles are in motion so they are changing the coordinates of the phase 
space as time goes by, but the density function in the equilibrium is invariant. 
Therefore the ergodic hypothesis is satisfied: the time average of any variable 
associated to the system coincides to its ensemble average. 

However, as the evolution in time of a static packing is not produced by 
itself, the distribution of the packing in the phase space can only be changed by 
driving the system (shaking, taping, compressing,...). It gives the opportunity 
of being unconcerned about whether solutions are steady in time or not or about 
which the time average of any quantity associated to a packing is. Then, follow- 
ing the Boltzmann approach, it is possible to compare different packings of a 
static and jammed granular system (microstatc) compatible to the same exter- 
nal constraints. An ensemble is considered as a set of copies of the packing that 
correspond to all different configurations for a given macrostate. Independently 
of how they had been arranged (i.e. independently of the dynamics during the 
driving process), they can be identified with a set of points in a suitable phase 
space. In the absence of further information there is not any a priori reason for 
favoring one of these microstates more than any other. Hence it should be nat- 
ural constructing the equilibrium function by assigning equal statistical weights 
to all the functions compatible with the requirements of the ensemble. 

Furthermore, when ideal gases are considered, the distribution function can 
be analyzed in the phase space of a single particle. But, in the end, this is 
possible because the Boltzmann factor, i.e. the statistical weight of the points 
of the phase space, depends on the Hamiltonian function, which is an addi- 
tive and completely separable function (i.e. J-{qi,q2, ■• • ,QN,Pi,P2, mmm ,Pn) = 
f(SLi,Pi))- It means that the possibility of using the phase space of one 
particle just lies in a mathematical feature, but it does not depend on how 
the mechanical interaction of the particles is. Similarly, approaching solutions 
of static granular packings compatible to a macroscopic constraint is formally 
analogous. Both the ensembles reported in section 2 above are also based on 
additive and separable functions: the total elastic energy and the total force 
moment of a microstate are equal to the sum of the separate contribution of 
each particle. Therefore, following the Boltzmann approach and using the r, T 
phase space, it is possible to analyze solutions as if the system were decoupled. 
Respectively, the angoricity and the average elastic energy play the role of the 
"temperature" , the total force-moment tensor or the total elastic potential en- 
ergy play the role of the "Hamiltonian function". They only depend on the 
"position" of the particles in the "phase space" (i.e. on coordinates or af^) 
but they do not depend on the "relative position" of the other particles (i.e. 
Y*ij — ■ or afj — 0y). The whole partition function is equal to the partition 
function of one particle to the power N and the "phase space" of a single parti- 
cle can be used. So that the number of dimensions of this phase space is small, 
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(K + l)K/2. It is very important to emphasize that this assertion is just based 
on mathematical arguments. As a result of the mechanical interaction between 
particles, they can be located in stress or force-moment levels, but for statistical 
purposes the system is decoupled. 

At this point it is possible to analyze how the nature of the phase space 
of a single particle is. The r, T phase space has two important features: 1) 
its accessible region (given by the integration limits of the coordinates and the 
volume element such as the points within this volume can be compatible with 
a physical solution of a given problem) and 2) the degeneracy of the points, 
or density of states (which comes from the capability of geometrical patterns, 
or configurations, to cover regions of the phase space). Both these issues are 
analyzed below together with an explanation for why using a change of variables 
can be interesting. 



4 Coordinates (p, q, uj) and (P, Q, Q) 

If the ensemble is based on the total elastic potential energy of the granular 
system, a natural workshop is provided by the components of the stress tensor 
Gij [9]. In the static case, for elastic continuum media the energy density is 
given by 

dE 1 o m 

where Sijki is the compliance tensor of the material. 

Within granular media the particles are in equilibrium under a finite number 
of forces. If they are elastic, they store elastic potential energy, which for a 
particle fc is given by 

E h = 2 J vk VijSijmn°Ln dVk = 9 j yk S ijmn a mn dV P — 7J ^2 ^ JC^J ' 

(2) 

where is the actual stress field within the volume of the particle Vp , Sfj mn 
is the compliance tensor of the material, cr& is the stress field averaged not only 
over the volume of the particle but also on the associated void (i.e., over the 
respective Voronoi cell V k ), and Sfj mn is the equivalent compliance tensor. f kl 
are the interacting forces 3 and K-q is the stiffness of the particles (i.e., the ratio 
between the applied force and the reduction of the diameter). 

In order to translate a set of forces into an equivalent average stress field the 

3 The final equality of the expression is only for linear force-displacement relationships, 
e.g., elastic disks. 
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average stress tensor [22] can be used: 



=fc f k j-t/k L „kl {kl 



(3) 



where r is the position of the point of application of the force / 



kl 



According to the possible arrangements, some set of forces can be equivalent 
to the same stress field. Therefore, the phase space determined by is a good 
workshop to compare microstates, and according to the value of the forces it is 
possible to calculate the elastic potential energy. This is the r phase space. 

Nevertheless, the nature of these media makes it preferable to use other 
stress state variables. For the case of 2D media, considering the invariants of 
the second-order stress tensor aij , the variables (p, q, w), (aj, an, oS), or (er, m, ui) 
are more suitable. 

The coordinates p, q, oj are based on Mohr's circle and are defined as 
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(4) 
(5) 
(6) 



p and q 4 (or q' = q/p), the centrum and radius of Mohr's circle, relate to the 
invariants I\ = an and I2 = tija X id V j of the Cauchy stress tensor (being 
the Levi-Civita symbol); therefore, they arc also independent of the reference 
system. However, u> depends precisely on the orientation of the principal direc- 
tions of the tensor. Other alternatives arc to use the principal stresses a\ and 
an (also invariants) together with lo, or the minor principal stress a and the 
anisotropy factor m (the ratio between both stresses m = ai/au). 

Using these coordinates rather than the components of the Cauchy stress 
tensor is interesting because some physical restrictions and some magnitudes 
are conveniently expressed as functions of them; for instance, the elastic po- 
tential energy of an isotropic continuum medium does not depend on u, but it 
often does when the medium is not isotropic (e.g., a particulate medium). 
In this paper, (p, q' , w) are used for the stress phase space. In short, they express 
the stress level, the relative deviatoric component of the stress state, and the 
orientation, respectively. 



4 This notation of stress as p, q is widely used in triaxial tests performed in soil mechanics. 
Obviously they are not the position and momentum in the T phase space. 
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The volume of the accessible region of the r phase space of one particle V T 
can be expressed as 

j.cxd poo ry/OasntCyy fOO rl p2lT 

V T = / / / da xx da yy da xy = / / / 2p 2 q' dpdq' dui . (7) 

JO JO J -y/o-nxVyy JO JO Jo 

The integration limits 5 are those which guarantee that, for any state, the re- 
spective principal stresses are positive. It means that the stress field is compres- 
sive in any spatial direction (as it is expected for cohesionless materials). The 
principal stresses are the eigenvalues of the Cauchy stress tensor and, accord- 
ing to the characteristic equation, both of them are positive when a 2 y < a xx a yy . 

On the other hand, the force-moment approach is based on the force-moment 
tensor E^ = a^V k = ^ r\ l fj l . For a system of N particles, the total Ey is 
given by the sum of the force- moment tensor of every particle E,y = E^. 
Therefore, the coordinates of the phase space can be directly E^ . It is being 
referred herein to this space as the T phase space. 

However, according to the invariants of Ey, it is also possible to define the co- 
ordinates P, Q, n as P = A/2, Q = \J (1/4)1? - I 2 , and = (1 /2) arctan s 

being I\ = E,i and I2 = eijT, X iT, y j the invariants of Ey. Q' is again the relative 
deviatoric component, Q' = Q/P, and the Levi-Civita symbol. 

In this paper, (P,Q',n) are proposed for the force-moment phase space. 
The integration limits and the volume clement are similar to those used above: 
P 6 [0, 00), Q' G [0, 1], n e [0, 2tt) and dV T = 2P 2 Q'dPdQ'dn. 



5 Canonical ensembles 

Consider the case in which the average external stress field is kept constant, al- 
though the system is driven many times to make it reach the equilibrium. The 
system is supposed to be in contact with an external pressure bath. 

Following the approach of elastic potential energy, the expected statistical 
weight of a configuration in equilibrium would be governed by a Boltzmann 
factor which depends on quadratic functions of the stresses as 

E c,ip2, q >2,u,) 

e 3 , ( 8 ) 



5 Actually, u> should be evaluated within the interval [0, 7r), because the stress states of ui 
and lu + 7r are the same. Nevertheless, this doubled evaluation results in a more intuitive 
representation. Note that, in this way, a determined system of forces, equivalent to a stress 
field, would be represented by two symmetrical points in the phase space. 



where f) relates to the external pressure and the expected equivalent stiffness of 
the system, and the energy depends on the configuration and on the stresses. 
This approach does not mean that the energy is constant during the driving 
process, but that, once the kinetic energy has been dissipated and the packing 
is again in mechanical equilibrium with the external pressure, the elastic energy 
of the ensemble is the same that it was before the driving. As the energies 
can be expressed as E a = p 2 ■ M a (g' 2 ,w) ( see [9], although the function was 
expressed in terms of other variables E a = Ei(K D ,a 2 )M a ( mjU> )), the probability 
of the configurations actually depends on the functions Mj^^j . The energy 
increases with the stress level and with the anisotropy factor, so the statistical 
weight decays. After imposing the restrictions of this approach, the following 
important conclusion is obtained: the probability of the configurations depends 
strongly on the volume of the respective allowed regions. This consideration 
can be useful for other approaches. 

On the other hand, following the force-moment approach, the Boltzmann 
factor would depend on the force-moment tensor [10, 11, 17, 18] as 

e -X«E W) ( 9 ) 

where Xij 1S the angoricity, a second-order tensor. Using a reference system 
x, y according to the principal directions of the angoricity, i.e., those in which 
Xxy = 0, the statistical weight does not depend on £ xy . 

The physical meaning of the angoricity can be deduced by obtaining average 
values. For instance, by integrating over the volume of the T phase space of a 
single particle (analogous of Eq. 7 for Sy variables, instead of it is obtained 
that the expected value of T, xx is given by: 

y f°° yi 3 / 2 p-£ Ix x**^y q 1 

XX N ~ J- S V V ^WF_ " 2 Xxx (10) 
The 3/2 and 1/2 powers on T, xx come from the integration limits of Ti xy (£ xy £ 

{~ y^xx^yyi + y/^xx^yy)) ■ 

Acting in the same way, < Y, yy >= (3/2)(l/x TO ) and < Y> xy >= 0. Using 
£jj = V<Tij and aij = r,j it is shown that Xij i s actually 3/2 of the inverse of 
the average external force- moment tensor VYy per particle, i.e., 

Xij = 2y{nj) , (if) 

where V is the volume of the whole system. 

However, the exponent of the Boltzmann factor can also be expressed as 

Xij^ij = Pxp [1 + Q' cos 2n XQ 'n] , (12) 



6 The reference system is established according to the orientation of the configuration, so 
that id refers here to the angle of the principal stresses with respect to the configuration. 
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where \p = (Xxx + Xyy) an d XQ'n = (Xyy - Xxx)/(Xyy + Xxx) & [0, 1). The 
advantage of writing the Boltzmann factor in this way is that, in the statistical 
weight of each point of phase space, it separates the contribution due to the 
force-stress level P from the contribution due to the relative deviatoric compo- 
nent Q' and the orientation SI, depending on the external control parameters 
Xp and XQ'n- 

As some physical restrictions are independent of the stress level, it is inter- 
esting to consider the distribution of points in Q' , SI planes, for any value of P or 
for all of them. The integration of the points over all the values of P and its rep- 
resentation over a polar Q 1 , SI plane is referred to as the comparison plane (i.e. 
a plot which displays Q' sin fi as ordinate plotted against Q' cos SI as abscissa) 7 . 

If every point in the T phase space of a single particle is supposed to have 
the same multiplicity, the partition function (when negative stresses are not al- 
lowed) is given by 

POO pi /"27T 

Z t (1, X p,XQ'q)=2 / / e-^ p ^+^'^' cos2 ^P 2 Q'dPdQ'dn. (13) 
Jo Jo Jo 

Integrating over all the values of P leads to the expected density in the com- 
parison plane P(Q>,n)i which just depends on XQ',Q- 

For the case of isotropic compression, Xxx = Xyy = Xp/^, so that the ex- 
ponent of the Boltzmann factor just depends on P, i.e., Xij^ij = Pxp- This 
means that, if the granular system is in equilibrium with an external isotropic 
stress field, the probability of the particles in the ensemble will not depend on 
either Q' or SI, so that the expected statistical distribution of particles according 
to their respective values of P will be 



Pt(p) — 




(14) 



satisfying p T r P ^P 2 dP = 1.0. Note that the volume element is dVrp = P 2 dP 
and it affects the expected number of particles rip in each P ± AP/2 interval: 

^ = & e -* pP P 2 AP (15) 

Moreover, integrating over P, in the isotropic compression case every point in 
the comparison plane is supposed to be equally likely to be found 

PT(Q',n) = -, ( 16 ) 
satisfying ^ p T{Q , M) Q' dQ' dSl = 1.0. 



7 In a cfij phase space, the comparison plane plots the polar distribution of q' and u), 
integrated for all p values. 
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On the other hand, for anisotropic stress fields, the higher the anisotropy 
ratio, the stronger the dependence of the Boltzmann factor on Q' and f2. This 
is shown in Fig. 1, where it is plotted for different values of XQ'ti an d P- As 
fi is fixed according to the principal directions of the external stress field (via 
the angoricity) , the expected local stress fields are not isotropically oriented but 
according to some prevailing directions. 




Figure 1: Representation of the Boltzmann factor (using Eqs. 9 and 12) over 
two constant- P planes according to the force-moment approach. \p — 1-0- The 
higher the P, the smaller the value of the Boltzmann factor, whereas XQ'n se ^ s 
the dependence on Q' and Q. (a) P = 1.0 and XQ'Ct = 0.33, (b) P = 1.0 and 
XQ'n = 1.0, (c) P = 2.0 and X Q'U = 0.33, and (d) P = 2.0 and X Q>n = 1-0. If 
XQ'O i s 0.0 (meaning that the external stress is isotropic), every point within 
any constant-P plane is equally likely to be filled and uniform gray circles are 
expected. 
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6 Configurations. Allowed and forbidden re- 
gions 

6.1 The degeneracy of the points in the phase space 

Irrespective of which canonical ensemble is used, it is highly relevant that some 
points of the stress phase space can be filled by several configurations (and some- 
times with several orientations and several internal degrees of freedom) while 
other points are forbidden for some of the configurations. Consequently, the 
phase space must be characterized by a density of states function which comes 
from the sum of the allowed regions of all configurations in all possible orien- 
tations with respect to the principal directions of the external pressure. Then, 
Eq. 13 becomes 

/>00 />1 />2-7T 

Zr(l,XP,XQ>n) = / / g {QI) e^ni+x Q ,nQ'cos2n] 2p 2 Q , dpdQld ^ 
Jo Jo Jo 

(17) 

where <?(Q') is the total density of states function, which is supposed to be in- 
dependent of both fl and P. 

Since this function is actually equal to the sum of the respective contribution 
of every configuration g a (Q'), it is also possible to express the partition function 

as 

/>00 />1 />27T 

Z T (l, X p,XQ'n)=J2 / / / 9 am e XpP[1+x ^ Q ' cos2n] 2P 2 Q'dPdQ'dn. 
a Jo Jo Jo 

(18) 

The functions g a (Q r ) ar e equal to within the respective forbidden regions 
but can take a normalized value within allowed regions (because some configu- 
rations could have also internal degeneracy). Expressing the partition function 
in this way reflects exactly not only the stress state in which particles of every 
microstatc arc located but also the configuration to which they belong. Then, 
some expected properties, especially the associated packing ratio, can be di- 
rectly measured if the distribution function is separable in configurations. 



6.2 Physical reasons for allowed and forbidden regions 

Packing spheres is inherently a geometrical problem due to exclusion-volume 
effects, and how particles can locally be arranged affects the expected distri- 
bution of the disordered total arrangement. Two conceptual approaches for 
the study of jammed packings have emerged: the "ensemble" aproach and the 
"geometric-structure" approach [23]. The elastic energy approach [9] includes 
features of both approaches, by considering the mechanical stability of geomet- 
rically compatible local packings within an ensemble. This procedure is also 
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followed herein. 



A configuration a. refers [9] to a fixed number of particles arranged in a spe- 
cific geometrical pattern such that the cluster is in equilibrium with the stress 
field. In this way, some crystal-like configurations can be analyzed. However, 
there are other possibilities for the definition of the configuration. In this work, 
they are distinguished according to the number of forces. 

The number of degrees of freedom and their stability determine the allowed 
region of the phase space, which depends on the kind of interaction between the 
particles and on the geometric compatibilities. An important advantage of using 
the coordinates (P, Q',Cl) is that these physical restrictions are independent of 
the force-moment level (the ability of a configuration to fill a point only depends 
on Q' and f2), so that the comparison plane can also be used to analyze the size 
and shape of the respective allowed region. 

To illustrate this assertion, the simplest case of a 2D frictionless monodis- 
perse system is analyzed. In this medium, only five configurations are possible, 
defined respectively by the number of forces L: 6, 5, 4, 3, and 2. The forces f kl 
acting on the particle are ordered according to the angle 9 kl with respect to a 
fixed reference system associated with the configuration. There are three kind 
of restrictions: 

• Restrictions due to the kind of interaction 

In cohesionless materials, all forces must be positive: 



As there is no friction, all forces are central, and there arc no moments. 

• Restrictions due to geometrical incompatibilities 

For a 2D monodisperse system the maximum number of forces is six, and 
there are also some geometrical restrictions: 



(19) 



• Restrictions due to mechanical equilibrium 
In 2D there are just two restrictions: 




(20) 



Vfc,Z6» 



,kl 







> tt/3 



(21) 



and one equation 




(22) 



with e k0 = e kL - 2tt. 
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As a result, every configuration has an internal number of degrees of freedom 
with bounded values of the respective coordinates. All these restrictions are in- 
dependent of the force scale (if the forces are increased by K times, restrictions 
given by Eqs. 19, 20, 21, and 22 are still satisfied), and in any constant- P plane 
the allowed regions extend over the same Q',Sl values. So, it is not necessary 
to analyze the whole volume of the phase space, but just the comparison plane. 
Irrespective of the configuration, if the force-moment state approach is followed, 
every allowed point has the same statistical weight, whereas following the clastic 
energy approach it depends on the value of the forces. 

Then, the respective g a (Q') of a configuration is equal to the normalized inte- 
gration over all possible orientations of the geometrical patterns with respect to 
the principal stresses (fixed by the principal directions of the angoricity) . These 
functions include the features of the particles (size distribution, shape, rough- 
ness, strength, etc.) because this determines which configurations are possible 
and which regions of the phase space are covered. 



7 Probability of configurations and transitions 

Although not only the statistical weight but also the density of states affect the 
expected statistical distribution, probably particles are not totally free to move 
between configurations, and not all transitions are likely to occur during any 
cyclic process. As a result, the expected state after a protocol could relate to a 
Markov process, determined by the allowed regions. This could explain why, for 
instance, in 2D monodisperse media, protocols often cause most of the particles 
to arrange according to the six-force configuration, although other configura- 
tions may also be expected (because they also have considerable allowed areas) 

So, consider a particle in equilibrium with some forces (i.e., belonging to 
a configuration) in a system of N particles. If the external conditions of the 
whole system (i.e., the angoricity) are modified, another solution of the elastic 
problem with a different arrangement could be achieved. During this process, 
the forces on the particle can change, and it is possible to draw its path in the 
phase space. If this path remains within the allowed region, the particle can be 
kept in the same configuration. The particle could also change into other con- 
figurations without breaking if they overlap. Finally, if the stress path reaches 
a region where the configuration is not possible, it will break locally. Once this 
happens, the particle could fall into any configuration, according to the proba- 
bility of each. 

As a consequence, the stochastic evolution of the system in a determined 
driving process could be supposed to be a Markov chain in which the proba- 
bility of a particle belonging to a configuration a in the (N + l)th step P& +1 
depends on the distribution in the current step N and on the transition proba- 
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bilities p a p . These could be associated with the features of the allowed regions 
of the configurations. If the probability of changing the stress stress is pi, the 
probability of breaking is p2 a , the probability of changing to overlapping con- 
figurations is P3ap, and the probability of falling into a configuration once the 
equilibrium has been broken is po a , it follows that 

P^=p a0 P», (23) 

where 

Paa = (1 - Pi) +Pl ■ (1 - P2a) • (1 - P3ap) + Pi ' P2a ' POa (24) 

and 

Pap = Pi ■ (1 - P20) ■ PZafS + Pi ■ P2/3 ' POa- (25) 

The probability of changing the stress state P\ would depend on the cyclic mod- 
ification of the external control parameters, and this could explain the reported 
strong dependency on the protocol. The probability of falling into a configura- 
tion P a Q could relate to the size of the allowed region on the comparison plane, 
whereas the probability of breaking P2 a could relate to the shape of the region. 

8 Qualitative analysis of a 2D monodisperse sys- 
tem 

8.1 Configurations and allowed regions 

The possible equilibrium configurations of a 2D monodisperse granular system 
have been analyzed. This was done by numerically generating combinations of 
forces and angles which maintain the equilibrium and which correspond to the 
same value of P (P = 1.0). Then, points were plotted on the comparison plane. 

The surveying of combinations was done by fixing the angles of each of the 
analyzed configurations and then taking possible values of some of the forces 
within the bounded interval. The number of variables depends on the num- 
ber of degrees of freedom of the configuration. The domain of the variables 
was divided into a uniform grid, according to a fixed interval, and the other 
forces were obtained by applying the restrictions of force balance. The result- 
ing groups which included negative forces were automatically discarded. The 
value of P = 1.0 was imposing just by adding one restriction more. This is not 
necessary but it saves computation time because it avoids considering groups 
of forces which are actually scalar transformation each other (they would be 
plotted at the same point in a polar representation of (Q' , ft)). The interval of 
the grid was taken with values of different orders of magnitude, and although 
the number of obtained points was obviously increased with smaller intervals, 
it was observed that all the points were always located within the same region. 
Moreover, when the elastic energy approach was set up [9] the boundaries of 
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the six-force configuration were analytically obtained (although other variables 
were used in that paper) and that procedure gave the same result. Some con- 
figurations analyzed in this way are presented. 

The six-force configuration has five degrees of freedom (four forces and one 
angle) because inequalities of angles (Eq. 21) become equalities. So, once the 
angle is fixed with respect to a reference system, it is possible to analyze which 
combinations of positive forces keeping the equilibrium are possible. Plotting all 
these points over the comparison plane reveals the allowed area of this configu- 
ration (Fig. 2). It is a key factor that there are three orientations with respect 
to the principal directions of the equivalent stress field in which Q' can take any 
value between 0.0 and 1.0. This is absolutely the most stable configuration. 

The five-force configuration has eight degrees of freedom (three forces and 
five angles), although angles are strongly bounded due to the inequalities. It 
is not clear if every combination of angles can be considered as an independent 
configuration, because the multiplicity of some (Q', O) points would be too high. 
In any case, some representative cases are plotted in Figs. 3, 4, 5, and 6. It is 
important to realize that there are only points with high values of Q' for a few ge- 
ometrical patterns. Indeed, the highest value of Q' = 1.0 is obtained precisely in 
the quasihexagonal arrangement (which coincides with a six-force configuration 
in which one force is null). Moreover, the shape of some areas makes it reason- 
able to consider that the probability of breaking this configurations is very high. 

The four-force configuration has six degrees of freedom (two forces and four 
angles). Angles are somewhat less bound than in the five-force case. Some 
combinations can cover high values of Q, but only over a narrow area (indeed, 
only some lines), so that this configuration is in general quite unstable. Again, 
it is not clear whether the angles can change without breaking, which would 
extend the stability area of the configuration. This idea was used in a previous 
paper [9], where it was supposed that a transition between rhombic configu- 
rations happens until the system reaches a totally forbidden area. The result 
would be a four-pointed star area, covering some high values of Q. 



8.2 Numerical simulations 

Two different numerical molecular dynamics simulations were performed to 
qualitatively analyze the distribution in the r, T phase space. The purpose of 
these simulations was checking whether, at different stages of the compression 
process, the distribution of the points on the comparison plane is affected by the 
allowed regions of configurations or not. However the macroscopic constraints 
of the ensembles explained above are not totally imposed, because it is much 
more easy to perform simulations in which the displacement of the walls or the 
external pressure are controlled, rather than the force-moment or the final elas- 
tic energy of the packing. Nevertheless, as it is not being measured neither the 
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entropy nor the ensemble average of any quantity, these simulations are valid 
for the objectives of this paper. In addition, after a cyclic compression schedule, 
if the variation of the volume of the sample is small, then the constraints of the 
force-moment ensemble are closely matched. 

The LAMMPS code, a parallel particle simulator developed at Sandia Na- 
tional Laboratory [24] was used (including the GRANULAR package). The 
systems consisted of 900 (30 x 30) and 4, 900 (70 x 70) frictionlcss disks 8 , 
which interacted via a viscous and linear (Hookean) contact law (i.e. F kl = 
(Kud — , y 1 £S)n kl , being 5 the overlapping between particles k and I). A normal 
damping coefficient 7 ~ 1\[2\J Kn/m was used in order to dissipate the kinetic 
energy and compare static solutions. Particles were initially arranged according 
to a regular square lattice. Some irregular walls were also generated, two of 
which were fixed while the others were subject to controlled motion or applied 
force. The compression schedule is described by a first isotropic compression 
stage and then some anisotropic compression cycles. The sole purpose of the 
irregularity of the walls was to avoid the whole crystallization of the sample. 
These walls were generated by locating particles in consecutive contact and in 
such a way that they separate from the average position of the wall according 
to a normal distribution. Both the duration of the isotropic compression stage 
and the period of the cycles of compression At were chosen in such a way that 
At/t c ~ 10 4 , being t c the characteristic time t c = ym/ Kn. In the case of 900 
particles, an external pressure of p/Kn ~ 10~ 3 was applied (p is given in units 
of force per unit of length) whereas in the sample of 4, 900 particles a displace- 
ment of the walls of AL ~ 0.07Lo was imposed (being Lq ~ 70 x D and D the 
diameter). 

Fig. 7 and Fig. 9 show the aspect of the granular packings, and Fig. 8 and 
Fig. 10 show how the respective population density on the comparison plane 
of the two simulations was. In the first case, this measure was taken after 4 
anisotropic compression cycles whereas in the second case it was taken at the 
end of the isotropic compression stage. The first conclusion is that, although it 
is an isotropic compression state, the points are not uniformly distributed within 
the circle as expected for an isotropic compression state, providing evidence that 
the Boltzmann factor is not the only parameter which determines it (at least if 
the force-moment approach is followed, because results show that, the lower the 
value of Q', the higher the density of points). Therefore, the density of states 
g(Q') strongly determines the expected distribution of particles in force- moment 
levels. 

The number of points obtained in these simulations is not enough to ensure 
thermodynamic limit conditions. Therefore, obtaining histograms (which would 
be necessary to know how the g a (Q r ) functions arc and how the distribution de- 

8 Spheres were actually used, instead of disks, in a 2D system. Nevertheless, some precau- 
tions were taken in order to assimilate this system to a real bidimcnsional system. 



17 



pends on the protocol) it is not still possible. Nevertheless, this plot shows how 
the allowed regions have a strong influence on the statistical distribution. It is 
quite interesting that many of the particles are precisely located at the bound- 
aries of the allowed and forbidden regions of the main six- force arrangement. 



9 Conclusions 

It is possible to compare different packings of a granular system compatible to 
the same external constraints, in a suitable phase space which considers the 
equivalent stress field a\j or the force-moment tensor of the particles S*j- . The 
accessible region of the phase space depends on the nature of the granular sys- 
tem. 

According to the scale invariance of the mechanical restrictions, for 2D sys- 
tems it is not necessary to analyze the whole one-particle phase space but just 
the P distribution and the area in the comparison plane [polar representation 
of (Q',fi) for all P}. 

Due to the nature of the equilibrium states in particles, some regions are 
allowed whereas others are forbidden, meaning that some points of the phase 
space are degenerate. It can be expressed as a density of states function. 

The size and shape of the allowed regions of every configuration depend 
on the degrees of freedom and on the geometrical incompatibilities, mechanical 
constraints, and the kind of interaction between the particles. 

The statistical weight of the points in the phase space is fixed by the en- 
semble. Nevertheless, the final statistical distribution is strongly affected by 
the density of states, what is explained according to the allowed and forbidden 
regions of configurations. The intrinsic properties of the medium, i.e., the fea- 
tures of the particles (size distribution, shape, roughness, etc.), determine these 
functions, whereas the external stress state fixes the angoricity or the elastic 
energy (the parameter controlling the statistical distribution). More research is 
needed to analyze gmn and how it can be separated according to the respective 
contributions of the different configurations g a (Q')- 

As granular systems do not tend to equilibrium by themselves, it is necessary 
to drive them. The kind of modifications applied determines the probability for a 
particle to modify its stress state and thereby transition between configurations. 
These probabilities could be related to the shape and size of the allowed regions. 
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Figure 2: Allowed area of a six- force configuration. The only possible combi- 
nation of angles corresponds to a regular hexagon. This is the configuration 
which covers the largest area in the constant-P plane. It is about 55% of the 
expected maximum volume (the circle, which corresponds to an ideal, infinite- 
strength continuum medium). The highest value of Q which can be resisted in 
six directions is P, and some wide areas of high values of Q are also covered 
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Figure 3: Allowed area of a five-force configuration with pentagonal orientation 
of forces. The area is quite isotropic, and its size is about 32% of the expected 
maximum volume. The highest value of Q which can be resisted in five directions 
is less than 0.6P 
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Figure 4: Allowed area of a five-force configuration with a slightly modified 
pentagonal orientation of forces. The area is slightly more anisotropic than in 
the regular pentagonal case. The size is about 31% of the expected maximum 
volume, and the maximum value of Q is higher than the regular case in just one 
direction 
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Figure 5: Allowed area of a mixed hexagonal-square configuration. This five- 
force configuration is usually found in the grain boundary of crystalline regions. 
The size is about 24% of the expected maximum volume, and the highest value 
of Q is again P, being possible in only one direction 
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Figure 6: Allowed area of a quasihexagonal five-force configuration, being equiv- 
alent to a six-force configuration in which one force is zero. The size is about 
19% of the expected maximum volume, and the highest value of Q is P. How- 
ever, it is possible in just two directions. As the shape is narrow and strongly 
anisotropic, this configuration is supposed to be quite unstable 
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Figure 7: Arrangement of 900 elastic disks after the isotropic compression stage 
and some anisotropic compression cycles. The sample was initially ordered 
according to a square lattice, but the boundaries were irregular. 
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(Q/P).cos(Omega) 

Figure 8: Distribution of points in the comparison plane of the simulation of 
isotropic compression of 900 disks after the isotropic compression stage and 4 
anisotropic compression cycles. The number of particles is far from the ther- 
modynamic limit conditions and most of the points are located within the re- 
spective allowed region of the corresponding six-force configuration associated 
to the main crystalline domain (Fig. 7). 
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Figure 9: Arrangement of 4, 900 elastic disks after the isotropic compression 
stage. The sample was initially ordered according to a square lattice, but the 
boundaries were irregular. In this case, two main crystalline domains (grains) 
are found in the arrangement. 
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Figure 10: Distribution of points in the comparison plane of the 4, 900-disk 
compression simulation. It is still clear that the respective allowed regions of 
the two grains affect the statistical distribution of the points: two overlapping 
six-pointed stars can be discerned. 
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